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Abstract 

Guided by the Jacobi’s work published the year before his death about the 
rotation of a rigid body, the behaviour of the rotation matrix describing 
the dynamics of the free rigid body is studied. To illustrate this dynamics 
one draws on a unit sphere the trace of the three unit vectors, in the body 
system along the principal directions of inertia. A minimal set of properties 
of Jacobi’s elliptic functions are used, those which allow to compute with the 
necessary precision the dynamics of the rigid body without torques, the so 
called Euler’s top. Emphasis is on the paper published by Jacobi in 1850 on 
the explicit expression for the components of the rotation matrix. The tool 
used to compute the trajectories to be drawn are the Jacobi’s Fourier series 
for theta and eta functions with extremely fast convergence. The Jacobi’s 
sn, cn and dn functions, which are better known, are used also as ratios of 
theta functions which permit quick and accurate computation. Finally the 
main periodic part of the herpolhode curve was computed and graphically 
represented. 
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1 Introduction 


The history of the dynamics of the free asymmetric rigid body started with 
Euler. In the Whittaker’s book of mechanics one finds several references to 
the original papers ra- 

At present time, the dominant source for a detailed and profound treat¬ 
ment of the dynamics of the rigid body is the Klein and Sommerfeld treatise, 
recently translated into English [8]. 

The purpose to contributing to this venerable subject is the drawing of 
the predicted motion of the body represented by the trace of the three unit 
vectors, forming the rows of the rotation matrix, portraying on the unit 
sphere close trajectories representing the main periodic behaviour of the rigid 
body. Actually, the third row coincides with the representation of the angular 
momentum vector as is seen from the body system in the intersection of the 
sphere of angular momentum, and the ellipsoid of constant energy. This 
trajectory of the angular momentum vector has been graphically represented 
in different publications, see for example the Bender and Orzag book [T]. 

On the contrary, to my knowledge, the other two rows of the rotation 
matrix, have not received the pictorial representation of its motion, notwith¬ 
standing this motion for the second row has similar properties of symmetry 
that the angular momentum whereas the first row has different symmetries 
that the other two rows. In any case each of the three trajectories are periodic 
of the same period. The second and third row are symmetric with respect to 
the same two coordinate planes in the body system of principal moments of 
inertia. The first row, on the contrary, is skew symmetric with respect to the 
same planes. In my opinion this behaviour is an essential knowledge of the 
physics of the free rotation in space far from acting torques, which deserves 
to be known and taken in account. 

In order to present those facts and draw accurately as needed those curves, 
one uses a minimum of properties of elliptic functions, just those which are 
efficient and precise. Everything supported by the magistral Jacobi’s work. 

Our purpose is to dilute the embarrassment expressed by the historian of 
mathematics E. T. Bell [2] who writes: 

’’The rotation of the rigid body, for example, yields numerous elegant 
exercises in the elliptic theta functions; but few engineers who must busy 
themselves with rotation have time for elegant analysis.” 

After Jacobi, the excellent mathematician K. W. T. Weierstrass con¬ 
tributed in extraordinary form to the theory of the elliptic functions. In the 
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cited treatises both in Whittaker as in Klein and Sommerfeld, the properties 
of the new functions introduced by Weierstrass are used. However in this 
paper almost nothing of the Weierstrass functions are included. 

The basis to ignore the Weierstrass’ work, in this document, is justified 
by the following text extracted from the Whittaker and Watson treatise of 
Analysis [IB] which refers to the mathematics of the rigid body dynamics as: 

,: This result determines the mean precession about the invariable line in 
the motion of the rigid body relative to its center of gravity under forces 
whose resultant passes through its center of gravity. It is evident that, 
for purposes of computation, a result of this nature is preferable to the 
corresponding result in terms of Sigma-functions and Weierstrassian Zeta- 
functions, for the reason that the Theta-functions have a specially simple 
behaviour with respect to the real period-the period which is of importance 
in Applied Mathematics-and that the g-series are much better adapted for 
computation than the product by which the Sigma-function is most simply 
defined.” 

The previous text from two recognized experts in elliptic functions sup¬ 
ports the fact that for computing and also for drawing accurately our tools 
are the best. One can be more explicit: the theta functions used here are 
written in terms of Fourier series with n-coefhcients which are proportional 
to a number g, lower than 1, with an exponent n 2 . In such computations less 
than ten terms are necessary to obtain a precision of 15 places, for numbers 
q different from 1 in 1 over a million, one exceptional case, without gen¬ 
eral interest. One should use the theta functions because its extremely fast 
convergence. 

Richard Bellman confirms this point without any doubt in his brief book 
on theta functions [3]. 

The main objective of this paper is to represent by drawings the motion 
of a rigid body formed by particles which relative mutual distances do not 
change in time. The center of mass is assumed fixed at the origin of coor¬ 
dinates in the inertial system. The position of particle i of the body in the 
inertial system r, is known in terms of the entries of the rotation matrix, 
because the existence of a coordinate system fixed to the body which will 
be named the body system, with the same origin of coordinates as the in¬ 
ertial system, at the rest center of mass, where all the coordinates of the 
rigid body a* are constants of motion. The rotation matrix R transforms the 
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coordinates from the body system to the inertial system 


r* = Ra, , (for all i). 


( 1 ) 


In what follows we will be concentrated in the drawing of the components 
of the rotation matrix, which to maintain the constant distance between 
particles, requires its rows to be formed by the components of three unit 
vectors, mutually orthogonal, that form a right tern, namely, the x product 
of the two first vectors of R is equal to the third row of the same. The 
ortho-normality of the rows of the rotation matrix is represented in matrix 
notation as 

R t R = RR t = E , (2) 

where E is the unit matrix. The super-index T to the right of a vector 
or matrix denotes the transposed vector or matrix. Therefore one has the 
inverse matrix of the rotation matrix is equal to its transposed matrix. 

The derivative with respect to time is denoted with a point on the function 
to be derived. The derivative with respect to time of the equations in (2), 
uncovers the existence of the antisymmetric matrices defining the angular 
velocity 
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The first definition corresponds to the angular velocity in the body system 


<j) = 



(4) 


which gives the derivative with respect to time of the rotation matrix in the 
useful form 

R = Rw x . (5) 

The notation with the x product is used because the product by the right 
of the antisymmetric matrix with any vector is equal to the x product of the 
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vector of the matrix with the same vector 



( a, \ 

( ° 

—UJ 3 

W 2 \ 

( «1 

<J J X 

a - 2 = 

UJ 3 

0 

-Wi j 

0-2 


\a 3 J 

\ -<^2 

u i 

0 J 

\ a 3 


(6) 


The angular velocity in the inertial system f2 will be important in the last 
part of this work because one pretends to draw the curve, called herpolhode, 
portrayed by this vector when one projects it on the orthogonal plane to the 
angular momentum vector. 

From its definitions (3) the angular velocities vectors are related by equa¬ 
tion 

fix = Ru; x R 1 = (Rw)x , Q = Ru;. (7) 

where we have used the geometric theorem: the x product of rotated vectors 
is equal to the rotation of the x product of the vectors. 

To know the rotation matrix one starts using the conservation of the 
angular momentum vector. The angular momentum vector (with respect to 
the center of mass) is defined by the sum vector of the mass by the x product 
of the position and the velocity 

J = xri = ^mi(Ra ;: ) x (Raj 

i i 


= ^mjRaJ x (Rw x a*) = R^m^a* x (oj x a*), (8) 

i i 

where the positions were used in terms of the rotation matrix; the time 
derivative of the rotation matrix as a function of the angular velocity, and 
the geometric theorem: the x product of the rotated vectors is the rotation 
of the x product of the vectors. Next the double x product is used to obtain 


J = R£ 


rrij 
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(9) 


where one has introduced the unit matrix E to take out as a common factor 
the angular velocity vector. In this way one finds the angular momentum 
vector of the inertial system as the product of the rotation matrix, the inertia 
matrix of inertia in the body system, and the angular velocity of the same 
system 

J = RIw . (10) 
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The inertia matrix X is the following constant symmetric matrix, which is 
assumed diagonal since the body system could be selected in such a way. 


X = Y J m i 

i 
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Quantities Ij are positive and are called principal moments of inercia. The 
body system where the inertia matrix is diagonal is called the system of the 
principal moments of inertia. In what follows all the vectors of the body 
system have components in this system of coordinates. 


2 Computing the angular momentum in the 
body system 

The angular momentum vector is conserved. This is a general theorem of 
mechanics: if no external torques are present, the angular momentum vector 
is conserved. The inertial system of coordinates is selected so that the angular 
momentum vector is directed along the coordinate axis 3. The constant 
magnitude of this vector is J 


J T = J(0,0,1). 


( 12 ) 


According to equation (10) the components of the angular momentum in 
the body system are L =Xux The magnitude of this vector is the constant J 
as it does not change with the rotation, but its direction changes with time. 
One denotes by u the unit vector in the body system in the direction of this 
vector 

L = Xu; = Ju . (13) 


Hence the vector L is rotated into the vector J and the vector u is rotated 
into the vector k: 

/ 0 


J = RL , Ru = k = 


0 

Vi 


(14) 


As the rows of the rotation matrix are mutually orthogonal the second equa¬ 
tion in (14) implies that u T is the vector forming the third row of the rotation 
matrix. One proceeds to compute this vector. 
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The equation of motion for L is obtained from the derivative with respect 
to time of the previous equation 

0 = R(L + u; x L), u = u x u;, (15) 

where it has been written the time derivative of the matrix R in terms of the 
angular velocity. This is the Euler equation for the free rigid body. 

In other hand, according to (13), the angular velocity is written in terms 
of vector L as 

w=X” 1 L. (16) 

We find the equation of motion for the components of the angular momentum 
vector in the body system 

L = LxT 1 L. (17) 

Which has the constants of integration energy E and magnitude of the an¬ 
gular momentum vector 

L t X ] L = 2 E , L t L = J 2 . (18) 


We turn the attention to the unit vector u in the direction of the angular 
momentum, dividing the angular momentum vector by its magnitude. 

One writes the inverses of the principal moments of inertia and the con¬ 
stant 2 E/ J 2 in the form 


— H -\- Pcj , 


2 E 
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(19) 


where H and P are defined in order to lower the number of independent 
parameters 
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( 20 ) 


9 \T{ 'il ii hh hh hh 


e i + e 2 + e 3 = 3/2 , (21) 


which implies the parameters ej will be known in terms of one angle 7 

ei = cos 7 , e 2 = cos (7 — 27r/3), e 3 = cos(q + 27r/3). (22) 

The constants of motion that for the angular momentum vector were 
functions of five parameters E, J, Ii, J 2 , / 3 ; for the components Uj of the 
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unit vector in the direction of the angular momentum are written just in 
terms of the two parameters 7 and eo 

u i T u 2 + M 3 = 1, &\u\ + e 2 u^ + €- 3 M 3 = eo • (23) 

These equations suggest to utilize the sphero-conal coordinates aq and a 2 for 
the components of this vector 

Mi = cn(a 2 , fc 2 )cn(ai, £q), 

u 2 = dn(a 2 , k 2 )sa.(a 1 , Aq), 

u 3 = sn(a 2 , fc 2 )dn(aq, ki). (24) 

The sphero-conal coordinates are used at present time. According to the 
NIST Handbook of Mathematical Functions [IB], the parameter f3 in 29.18.2 
of this handbook corresponds in our notation to K{k\) + iK(k 2 ) — ia 2 , in 
terms of the function and the arguments which will be defined in this section. 
One prefers the form (24) since it is not necessary this far of paper to use 
complex variable. In the quantum case of the same Euler free rigid body 
the Schrodinger equation is separated in sphero-conal coordinates [12]. (See 
for example the R. Mendez-Fragoso and E. Ley-Koo ra review paper on 
quantum rotations.) The sphero-conal coordinates in the form (24) are also 
found in the Morse and Feshbach book HU, as conical coordinates. The 
stereoscopic view in 3 dimensions of some of the conical coordinate curves 
have been drawn in this reference. In our forward figure 1 some coordinate 
curves ( a 2 = constant) are drawn on the sphere as the trajectories followed 
by the angular momentum vector in the system of the principal moments of 
inertia. 

The Erst equation of (23) is satisfied identically with these coordinates if 
we use the properties of the Jacobi elliptic functions 

sn 2 (/3, k) + cn 2 (/3, k) — 1, & 2 sn 2 (/3, k) + dn 2 (/3, k) — 1 (25) 

and the relation of the constants ki 

k\ + k\ = 1 . (26) 

The second equation of (23) is satisfied identically with these coordinates 
if a 2 , ki , k 2 are constants chosen with the restrictions 

sn(a 2 , k 2 ) = J Cl 6 ° , cn(a 2 : k 2 ) = J C ° 63 , dn(a 2 ,k 2 ) = 

V e l — e 3 V e l — e 3 









(27) 


(gl ~ e 3)( e 2 ~ £o) 

(e 2 - e 3 )(ei - e 0 ) ’ 

which are redundant because the equations (25) and (26) are identically 
satisfied. Note that these constants are only functions of the parameters 7 of 
asymmetry and eo of energy. For the components u 1 , w 2 and M 3 of vector u, 
and of the third row of the rotation matrix R, it is preferable replace instead 
of (24) 

Ui = J C ° e 3 cn(oq, k x ) , 

\ e 1 - e 3 

u 2 = J C ° e 3 sn(ai, ki ), 

V e 2 - e 3 

dn(o; 1 , ki) . (28) 

Coordinate op is a function of time. To know its behaviour we must know 
the derivatives of these Jacobi functions. Owed to the quadratic relations 
(25), it is sufficient to know one of them, but one includes the three to 
remark its simplicity and similarity 

dsn(fi,k) 

-—-= cn(/3, k) dn(/3, k ), 

d k ^ = -sn(/3, k) dn(/3, /c), 

</(l ^ 1 = ~k 2 sn(f3, k) cn(/3, A:). (29) 

We might know now the equation of motion for the coordinate «i. Sub¬ 
stituting (13) in the Euler equation of motion, with the explicit form (28) for 
the components of u and the principal moments of inertia in the form (19). 
One finds the constant derivative 

di = PJ y f(ex — eo)(e 2 — e 3 ). (30) 

Note that besides the used parameters 7 and eo, the time appears without 
dimensions in the combination JPt. 

To draw the curve u, going over the unit sphere, we should learn to 
compute the Jacobi elliptic functions. As functions of the real variable op 
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they are periodic. The period of the vector u is 4 times the function K(k\). 
The function K(k) is a complete elliptic integral of first kind. To compute 
it we have the algorithm of the arithmetic-geometric mean [Ij allowing to 
calculate it with rapidity and is described as follows: Take the initial values 
x Q = 1 + k, y 0 = 1 - k. Iterate x n , y n ->• x n+i = (x n + y n )/ 2, y n+1 = y/x n y n , 
up to the desired precision the two averages are equal to a certain value 
x. Then K{k) = n/{2x). To draw the vector u one divides the period in 
equal parts and one draws the components of u calculated in all the points 
corresponding to the increments of the division of the period. 

We have several efficient algorithms to compute numerically the Jacobi 
functions |T]. In this paper one uses the Jacobi formula [7] as the ratio 
of Q(f3,k) and H(f3,k ) functions, which has the advantage of a very fast 
convergence, comparable in precision and fastness with other methods, and 
the bonus of being one of them, indispensable to the calculus of the other six 
entries of the rotation matrix R. 

The theta Jacobi functions require other function of the parameter k 
defined in terms of K{k) by 

q{k) = exp(— 7 iK (\/l — k 2 )/K(k ), (31) 


which is a positive number lower than one, which is shortened as q. The 
theta functions are computed then [7] as Fourier series of real period 2 K(k) 

0ftfc) = l + 2B-l)Yco.^, (32) 


and period 4 K(k) 

H(P, k) = 2f;<—l)” 9 <" +1/2|a senffi^±lM. (33) 

n =0 \ K ) 


These series should be computed up to the necessary precision when the next 

2 

term is negligible, because the fast convergency to zero of the factor q n . 
These theta functions are the algorithms to compute other elliptic Jacobi 
functions by means of the equations [7j 
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(35) 
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Figure 1: The sphere of angular momentum intersected by ellipsoids of con¬ 
stant energy. The curves follow the motion of the angular momentum vector 
on the sphere of angular momentum for different values of the energy. These 
curves are symmetric with respect to two of the coordinate planes. 
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Jacobi [7] provides the entries of the rotation matrix as ratios of these 
theta functions. He understood the computational advantage of their func¬ 
tions. In this work one has adopted some priority to the Jacobi’s elliptic 
functions sn, cn, dn whenever it is possible because they are better known 
functions and their properties. But the Jacobi’s theta functions are actually 
used as algorithms to compute the previous functions. Comparing the Ja¬ 
cobi’s expressions [7] for the third row of the rotation matrix we can observe 
the same components with the same functions. We are following Jacobi very 
closely. 

As an example of the use of this algorithm to draw on the angular mo¬ 
mentum sphere the intersections with ellipsoids of constant energy one has 
Fig. 1, in which the asymmetry parameter was selected as 7 = 7 t /5 and sev¬ 
eral energy parameters of energy e$. This is the trajectory followed by the 
unit vector along the third row of the rotation matrix on the sphere of unit 
radius. The observer of the picture has the spherical coordinates 27r/T, 27 t/ 7 
with respect to the principal inertia directions. 

3 The entries perpendicular to the body an¬ 
gular momentum in the rotation matrix 

In this section one analyses the paleography of the rotation matrix as a func¬ 
tion of time, of the motion without torques of a rigid body, as published by 
C. G. J. Jacobi [7J at 1850. Understood paleography as the study of an old 
text and its traduction or explanation in modern terms. The nine compo¬ 
nents of the rotation matrix were written by Jacobi by means of a set of 
theta functions, the same functions used in the previous section. The third 
row of the rotation matrix, one pointed before, was written in sphero-conal 
coordinates in terms of the H and 0 theta functions. Jacobi finds the other 
six components of the same matrix expressed also by using ratios of the same 
functions, although now with complex arguments. The corresponding entries 
in the two first rows of the rotation matrix involves three theta functions, 
repeated with different argument, and multiplied, by the real and the other 
by the imaginary part of a theta function with argument the complex vari¬ 
able an + za 2 - Remember an is proportional to time, and a 2 is a constant 
coordinate satisfying three equations in (27), but its numerical value should 
be determined now. 
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Figure 2: The trajectory followed by the vector forming the first row of 
the rotation matrix in a period. Perspective view and projection on the 
coordinate planes in the system of principal moments of inertia. The curve 
is antisymmetric with respect to the two coordinate planes crossing at the 
OZ axis of coordinates. 
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Since a 2 is an elliptic integral of first kind 


Ot-2 — 



dx 

\f{ i - - 1 - 1 ) f i - fc|^ 2 ) 


( 36 ) 


To attain a more precise value for ck 2 > one introduces the numerical value of 
this integral as approximated root and look for the root a 2 of the equation 


sn(a 2 , k 2 ) - \ — -- = 0 . (37) 

V ei - e 3 

Actually one finds two roots for ck 2 hi the interval [0, 2 K(k 2 )\. The smaller 
of the two, to be selected, correspond to the choice of positive values in (27). 

The 0 functions of complex argument are defined by the same Fourier 
series as before in (32) and (33). The modification due to the imaginary ar¬ 
gument changes the coefficient of the Fourier series, however the convergence 
remains very fast. 

Observing the two first rows of the rotation Jacobi’s matrix one deduces 
the convenience of combining the vector s of the first row as the real part, 
and vector t of the second row as minus the imaginary part of a complex 
vector, function of a complex variable. The two vectors are orthogonal to the 
third row, so the complex vector is too. The scalar product of the complex 
vector with itself should be zero as a consequence of the orthogonality of 
the real and imaginary parts and the same magnitude. We call null the 
vectors with square zero and note is not required the unit magnitude of the 
real and imaginary vectors to be null, it is sufficient orthogonality and same 
magnitude. 

It appears a common factor of the complex vector which is written in our 
notation as 


s — it 


H(ia 2 + K( fc 1 ))0(a 1 ) ^ 


Q(K(ki))H(a 1 + ia 2 ) \ 

—0(O)i7(cei + ia 2 + K{k\)) I , 
MH(K(ki))Q(ai + ia 2 ) ) 


(38) 


where one has suppressed the second argument of the theta functions because 
it is always the same: k\. One has changed also the sign of the first component 
because one assumes a Jacobi’s mistake to be verified in the sequel. 


14 



One notes that components 1 and 2 divided by the third allows to in¬ 
troduce the functions sn and cn of complex argument as the ratio of theta 
functions, 


Q(K(k 1 ))H{a 1 + ia 2 ) 
H(K(ki)) ©(«! + ia 2 ) 


sn («! + ia 2 , ki) 
sn(K(ki),k 1 ) 


sn («i + ia 2 , ki ), 


0(0) H{a i + ia 2 ) + K{k x )) 
H(K(ki)) 0(aq + ia 2 ) 


cn (aq + ia 2 , ki ) 
cn (0, ki) 


cn (aq + ia 2 , k ±), 


after this, appears a common factor given by the third component, then the 
complex vector can be written as 


s — it 


if(tf(fci))e(ai+ia 2 ) 

H{K(k l ) + ia 2 )Q{a 1 ) y 


sn(aq + ia 2 , Aq) \ 
•®cn(aq + ia 2 , ki) 

I 


(39) 


One notes the vector on the right is a null vector as a consequence of the first 
quadratic identity in (25), although the magnitudes of the real and imaginary 
parts of this vector are not equal to 1. 

One substitutes the equations discovered at page 575, equation 16.21 of 
the Abramowitz and Stegun handbook p] for functions sn and cn of complex 
argument, which have been written in function of the three components of 
the unit vector along the angular momentum. These formula were not found 
in other classical references and therefore were verified starting from the 
addition formulae for the Jacobi elliptic functions and the Jacobi equations 
for the same functions of imaginary argument which are reproduced in many 
references. 

u 2 + iu^Ui \ 

1 — 

—Ui + iu 2 u 3 
1 -u\ 

J 



H (K (ki))Q(ai + ia 2 ) 
H(K(ki) +ia 2 )0(ai) 


We verify in this form the print mistake in the first component of the com¬ 
plex vector which was corrected to be orthogonal both its real part and its 
imaginary part to the vector u, as it should be. 

To establish clearly these equations one rewrites them in terms of Euler 
angles as defined by Goldstein p]. Our matrix R corresponds to his A. 
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Figure 3: The trajectory followed by the vector forming the second row of 
the rotation matrix in a period. The curve is symmetric with respect to the 
same two coordinate planes that the body sistern angular momentum vector. 
Perspective view and projection on the coordinate planes in the system of 
principal moments of inertia. 
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The third row is the unit vector in the angular momentum direction which 
components are written in terms of two Euler angles 

u T = (ui,u 2 ,u 3 ) = (sin# sin0, sin# cos0, cos#). (41) 

Adding the first row to the second multiplied by the square root of minus 
one of the rotation matrix, it is possible to extract a common factor of the 
three complex components to find the null vector 

e*^(cos -0 + i sin 0 cos #, — sin 0 + i cos 0 cos 9, —i sin 9) . (42) 


It is interesting to observe that if the vector between parentheses is divided 
by sin 9 it appears again the vector one has found in the Jacobi’s publication 


/ 

V 


COS 0 

sin# 

sin0 

sin 9 


+ i sin 0 cot 9 

+ i cos 0 cot 9 

—i 


\ 

I 


( u 2 + iu 3 ui \ 
1 — u 3 

—u\ + iu 2 u 3 

l-«§ 

V -i 


( sn(«i + ia 2 , fci) \ 
—cn(o; 1 + ia 2 , k\) 

V -* / 


(43) 

For computations we should use the middle form because it has been 
written in terms of the components before computed of the unit vector u. 
What means an economy in computing machine time. 

The factor that should multiply this vector to obtain the two first rows 
of the rotation matrix as its real and imaginary parts of the vector (42) 
is e*^ 1 sin 9, where 0 1 is the main periodic part of the third Euler angle 0. 
According to Jacobi that factor is 


id, i U _ 2 = ^(A0fci))0(ai + id 2 ) 
V % H{K{k x ) + ia 2 )®{cL X ) 


cn(d 2 , k 2 ) 


©(O)O(ai + ia 2 ) 
0 (ai) 0 (zd 2 ) 


Actually, ploting the real vs. the imaginary parts of the expi0!, the points 
do not fall on a unit circle as it should be. The reason seams subtle to 
me. A graphical solution was to use the other a 2 root in the right hand 
side of this equation, which was denoted as a 2 which is equal to 2K(k 2 ) — 
a 2 . Modification of the right hand side is performed by the prescription 
16.33.4 in reference [1] who adds a constant angular velocity along the angular 
momentum vector with the same period that vector u 


J4> i. pr~2 = ^(^©(oq + m,) ( -umA 

V 3 HiKikJ+iaJGiaJ P \2K(k 1 )) 
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= cn(a 2 , k 2 ) 


0 (O)0(o!i + ia 2 ) 

0(ai)0(ia 2 ) 


exp 


—inai ^ 
2K(ki)) ' 


(44) 


This Jacobi’s equation allows to obtain the main periodic part of the 
third Euler angle which remains to know to complete the three Euler angles 
as a function of time. In fact one has the real quantities \Jl — «§, cn(a 2 , k 2 ), 
0(0), 0(ia 2 ), 0(«i), therefore taking the logarithm of (44) and subtracting 
its complex conjugated expression give us 


2 i(pi = 


0(ai + ioi 2 ) 

In ——- 

0 («i — ia 2 ) 


’’Kih) 


(45) 


as is found in references [14] and [9J, as a periodic part with the same period 
of the u vector. But Jacobi’s equation (44) contains more information since 
provides with economy of computation, the factor that is needed to compute 
the two missing rows of the rotation matrix. 

The difference of the two angles 0 and 0i is an angle with constant angular 
velocity. This was neglected by Jacobi due to its simplicity. One draws the 
figures of the unit vectors forming the Jacobi’s rotation matrix neglecting 
the same constant angular velocity around the constant angular momentum. 

In figures 2 and 3 are drawn the curves that follow those two rows of the 
rotation matrix on the unit sphere in a period AK(ki) of a±, for a typical 
particular case, computed from equation (44). 


4 Recovering the Jacobi’s expression for the 
components of the rotation matrix 

In this section one proofs the Jacobi’s equation (44) of the precedent section 
which could be used as an efficient algorithm to compute the drawing of the 
other curves of two rows of the rotation matrix. The full rotation matrix 
including the missing constant angular velocity rotation around the constant 
angular momentum vector is fully rcincorporated here. 

One presents a simplified mathematical proof to verify equation (44). As 
an alternative one could use equation (45) that seems to be more elegant and 
simple. However, as has been pointed out, it hides information, it implies 
first compute the angle, and then compute the trigonometric functions and 
replace in the rotation matrix. The Jacobi’s equation (44) give directly the 
trigonometric functions and permits ’’quick and accurate computation” [3]. 
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The derivative of the angle 0 is obtained by Landau and Lifshitz [9] from 
the angular velocity and angular momentum in terms of the Euler angles (the 
same convention for the Euler angles that used Goldstein | 6 j and is adopted 
in this work) 


OJiUi + UJ 2 U 2 


Un 


U U) 


U3LJ3 


U o 


2 E/J-ujJ/h _ J JP(e 0 -e 3 ) 

l-ul I3 1 - U3 


J_ 

h 


+ JPy/(ei - e 0 )((e 2 - e 3 ) 


eo — 63 


(ei - e 0 )((e 2 - e 3 ) 


1 — Un 


(46) 


The first term of the right hand side of this equation is a constant, with 
an integral lineal in time. The second summand is a periodic function of 
the variable aq with period 4iL(Aq), considered by Jacobi. The integral with 
respect to time of a periodic function is other function lineal in time plus a 
periodic function of the same period. The behaviour of 0 is then separated in 
two terms, one periodic function of the variable aq of period 4 K(k\) (the main 
peridic term) plus other function lineal in time, where the two contributions 
lineal in time are added. The last is a motion with constant angular velocity, 
with other period generally incommensurable to the previous one. This last 
with a constant angular velocity is separated and neglected by Jacobi due to 
its relative simplicity. In what follows this terms are not longer ignored. 

From the two terms in the right hand side of (46), the second which will 
be denoted dfo/dt, is convenient to express in terms of the coordinate aq 


<J0o _ eo — e 3 1 

**1 J(e 1: — *o)((ej - e 3 ) 1 - sn 2 (a 2 , k 2 )dn 2 (a u k t ) ' 


(47) 


The same derivative is now transformed replacing the constants in terms 
of elliptic functions of the imaginary value m 2 , which is accomplished by 
means of our equations (27), and Jacobi’s formulas for imaginary argument 

ra, la- 


d<t >0 .cn(ia 2 , k 1 )dn(ia 2 , ki) 1 

da± sn (iat 2 , ki) 1 — /cisn 2 (f« 2 , fci)sn 2 (aq, k\) 


(48) 


In what follows one suppressed the argument k\ since it is the same for all 
the functions. 


19 














It is customary to present the elliptic integral of third kind in the standard 
form [EJ 


ra i 

n(ai,ia 2 )= / da 
Jo 


ai A;jW(ia 2 )cn(ia 2 )dn(ia 2 )sn 2 (a!i) 


1 — /cfsn 2 (m 2 )sn 2 (ai) 


(49) 


which is close to the angle </> 0 (ai) because 


n(o!i, ia 2 ) 


cn(ia 2 )dn(ia 2 ) 

sn(za 2 ) 



da i 


1 

1 — /cfsn 2 (iQ; 2 )sn 2 (a;i) 


— aq 


(50) 


One finds so 

, . . . cn(fcK 2 )dn(ia 2 ) 

0o(«i) = «n(ai,za 2 ) + * -tt— 7 --aq , (51) 

sn(za 2 ) 

where the constant of integration is assumed to be 0 o (O) = 0. However 
instead of using the elliptic integral of third kind only the Jacobi’s Theta 
functions are used here. 

To obtain the expression of the factor used by Jacobi in his rotation 
matrix one finds a relation borrowed from the Whittaker and Watson book 
of analysis [15], p 518 


@'(aq-H ck 2 ) 0'(«i) B'(*o: 2 ) 
@(aq + ia 2 ) 0(aq) Q(ia 2 ) 


—fc 2 sn(aq)sn(ia; 2 )sn(aq + ia 2 ). (52) 


Which right hand side is written with the addition formula of Jacobi’s func¬ 
tion sn(aq + ia 2 ) used in (43) as 


0 '(aq + ia^) @'(aq) 0'(ia 2 ) 

@(aq + ia 2 ) 0(aq) 0(*a 2 ) 

—A;Jsn(m 2 )cn(ia 2 )dn(zQ; 2 )sn 2 (ai) 

1 — /cfsn 2 (aq)sn 2 (ia; 2 ) ~*~ 

A) 2 sn 2 (ai 2 , fc 2 )sn(aq, Aq)cn(aq, Aq)dn(aq, hi) 

1 — sn 2 (ai 2 , fc 2 )dn 2 (aq, hi) 

Integrating both sides of this equation from 0 to ay one finds 

@(O)0(aq + ia 2 ) ®'{ia 2 ) _ 
0(ai)0(za 2 ) 1 0(ia 2 ) 
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r a 1 — /c^sn(ia2)cn(ia 2 )dn(iQ;2)sn 2 (Q;i) \J 1 — sn 2 (a2, fc 2 )dn 2 (a:i, k\) 

Jo 1 — /ciSn 2 (o; 1 )sn 2 (iQ; 2 ) cn(a 2 , k 2 ) 

(54) 

where one replaces the angle 0o to find a version of the Jacobi equation (44) 
if one suppresses two of the three terms lineal in the time that are neglected 
by Jacobi 


0(O)0(ai + ia 2 ) _ . 
11 “* 


0i(o;i) + 


nai 


+ln 


\J 1 - sn 2 (a 2 , fc 2 )dn 2 (ai, h) 


2K(k 1 

where the lineal terms in time are no longer present as 

irai 


cn(a 2 , k 2 ) 


0 o(®i) — 0i (cq) — 


2K(ki 


+ ia.\ 


Q'(ia 2 ) cn(ia2)dn(ia 2 ) 


0(ia 2 ) 


sn(fa 2 ) 


(55) 


(56) 


This constant angular has been ignored in this work, as in Jacobi’s, taking 
in account only the main period of the rotation matrix, in which we have 
suppressed such constant angular velocity. One regards preferable this choice 
as discussed at once. 

That constant angular velocity has two terms, one of them depends on 
parameter H appearing in (20) and disappearing from discussion until reap¬ 
pears as a summand at the velocity of the angle of Euler 0 in (46): 


y- = JH + JPe 3 

h 


(57) 


One adds all the lineal terms in time to have the complete angle 0 minus 
its periodic 0 1; considered by Jacobi 


— 0i = JHt + tti 


e 0 


_ Q'( a 2; k 2 ) 7r(q 2 + A (k- 2 )) 

(ei — eo)(e 2 — e 3 ) ®( a 2 ik 2 ) 2K{ki)K(k 2 ) 


(58) 

The term proportional to H could have any value and give an undeter¬ 
mined character to any draw of the rotation matrix, as a consequence it is 
a vagary to pretend to draw it. Nevertheless it deserves some extra con¬ 
sideration. The first term and the first inside the parentheses should be 
reassembled as 

e o 


JHt T ol 1 


- E t 


\J(e 1 - eo)(e 2 - e 3 ) J 


(59) 
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In the remaining terms in that Euler angle we see the logarithmic deriva¬ 
tive of @(02) which is known in analysis as the Z(a 2). 


E 


(f) — (pi — —t-\-OL 1 


Q'{a 2 ,k 2 ) n(a 2 + K(k 2 )) 
0{a 2 ,k 2 ) + 2K(k 1 )K(k 2 ) 


(60) 


The logarithmic derivative Z(a 2j k 2 ) = Q'(a 2 ,k 2 )/Q(a 2 ,k 2 ) is computed 
as the ratio of Fourier series or it is used the algorithm of the arithmetic- 
geometric mean as is found in the literature pQ. 

As one has computed the nine entries of the main period of the rota¬ 
tion matrix, it is possible to draw on the unit sphere other three curves 
corresponding now to the three columns of it. They give the main periodic 
motion in the inertial system of coordinates of the rotating rigid body with 
the constant angular velocity around the constant angular momentum. The 
resulting curves are different to the previous result: one Ends one curve al¬ 
most plane symmetric with respect to two planes. The two other curves, are 
both symmetric with respect to one of the two coordinate planes. Trans¬ 
forming one into the other by reflection in the other two planes. The three 
are curves without crossing point. 


5 Drawing the herpolhode 

The equation of energy conservation when the components of the angular 
momentum and angular velocity vectors are given in the inertial system, im¬ 
ply that the component of the angular velocity along the angular momentum 
is constant 

2E/J = k T Cl. (61) 

The components of f2 perpendicular to k move on the plane perpendicular 
to k, called the invariable plane; its trajectory on the invariable plane follow 
the curve called herpolhode. The vector Cl minus the component of Cl along 
the angular momentum k is written with the double x product as 

(k x O) x k . (62) 

This vector in the body system is 

(uxu;)xu = uxu, (63) 

where equation (15) was utilized. 
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As a consequence the vector describing the herpolhode is 


R(ii x u). 


( 64 ) 


When studying the sphero-conal coordinates (24) one finds the tangent 
vectors to the coordinate lines 


e i = 


<9u 
da i 


( — cn(a 2 , &2)sn(cKi, &i)dn(a;i, k\) 
dn(a; 2 , A; 2 )cn(cKi, &i)dn(a;i, ki) 

\ — sn(a 2 , k 2 )kfsn(ai, fci)cn(ai, k\) 


(65) 


and 

-sn(ce 2 , k 2 ) dn(a 2 , k 2 ) cn(ai, ki) \ 

—/c|sn(a; 2 , k 2 )cii(a 2 , k 2 )sn(an, ki) . (66) 

cii(a 2 ,k 2 )dii(a 2 ,k 2 )dn(ai,ki) ) 

These vectors are perpendicular to vector u, and mutually perpendicular. 
The two vectors are of the same magnitude. Hence one has 



ei x u = —e 2 . 


(67) 


Note vectors u and ei are different by the constant factor di, therefore 
vectors u x u and —e 2 are different by the same constant factor. Besides, 
because ck 2 is a constant coordinate, one makes the substitution of the con¬ 
stants (27) in the vector e 2 . one deduces that the entries of the herpolhode 
Oi and fi 2 are 

( Ql \ 

H 2 = -diRe 2 , (68) 

V o ) 


Eliminating a common constant factor of the components of this vector one 
deals 

( K e i - e o)(e2 - e 3 ) , . . \ 

- cn(cKi, ki) 


R 


62 — 6 q 


V 


'(ei ~ e 3 )(e 2 - eg) 

6l — 6o 

/(e 0 - e 3 )(e 2 - e 3 ) 


sn(ai, k\ 


dn(o; 1 , k\ 


(69) 


) 


62 — 60 

We can trust to the computer the drawing and computation of the her¬ 
polhode until the rotation with constant angular velocity around the angular 
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Figure 4: The herpolhode without the term of constant angular velocity of 
the rotation matrix. It is tangent to two concentric circles. 


momentum vector. In Fig. 4 one finds this curve where one notes the few 
known properties of the herpolhode. It is symmetric with respect to one of 
the coordinate axis. We find the curve between two concentric circles. It 
is tangent to them at points on the coordinate axis of symmetry. As time 
increases, the tangent vector to the herpolhode rotates always in the same 
direction, until the value 2tt is reached in a complete period corresponding 
to 2 K(ki). Some discrepancy is apparent when comparing with pictures of 
the herpolhode in the literature [5j, but the motivation to this discrepancy 
lays in the fact that other authors do not suppress the action of the constant 
angular velocity that we ignore. This is an useful addition to the previous 
work in the literature. 
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